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Abstract 

This paper considers the extreme type-II Ginzburg-Landau equations, a nonlinear 
PDE model for describing the states of a wide range of superconductors. Based 
on properties of the Jacobian operator and an AMG strategy, a preconditioned 
Newton-Krylov method is constructed. After a finite-volume-type discretization, 
numerical experiments are done for representative two- and three-dimensional 
domains. Strong numerical evidence is provided that the number of Krylov 
iterations is independent of the dimension n of the solution space, yielding an 
overall solver complexity of O(n). 
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1. Introduction 

The nonlinear Schrodinger equation is used in many areas of science and 
technology and describes, for example, the propagation of solutions in fiber 
optics |T] and Bose-Einstein condensates in ultra-cold traps [2]. Prototypical 
for this type of models is the Ginzburg-Landau problem, widely used to study 
the state of both low- and high-temperature superconductors. Due to its highly 
nonlinear nature, the involved energy landscape, and the strong dependence of 
solutions on external conditions, numerical simulations of the Ginzburg-Landau 
model have become an essential tool for providing better insight into properties 
of superconductivity phenomena. 

The Ginzburg-Landau model has attracted wide interest since its inception 
in the 1950s. In particular, the work on the linearization by Abrikosov of the 
model around the upper critical field is widely known [SJ. The mathematical 
foundations for the equilibrium Ginzburg-Landau models are well developed 
^0 [5] and a framework for finite element and finite volume discretizations was 
provided [6j. Different types of discretizations and numerical approximations of 
the Ginzburg-Landau models have been developed since, all of which subject to 
numerical simulations. 

Throughout the physics literature, several methods for solving the Ginzburg- 
Landau equations are described. Used most prominently is a Gauss-Seidel-type 
iterative scheme [3 [8] that is readily implemented, yet fails to converge for 
systems with physically unstable vortex configurations. Furthermore, it only 
yields linear convergence close to a solution. In computational physics in general, 
the use of Newton-Krylov methods and nonlinear multigrid schemes such as 
FAS is widespread 0. Also, preconditioned Newton-Krylov methods are already 
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Figure 1: Typical solution ijr. f2 H> C (displayed as |^| 2 , argVO of the extreme- 
type-II Ginzburg-Landau equations, here for a flat triangular domain with 
circumradius 5 and the magnetic vector potential A(x, y) = (— y/2, x/2) T . Three 
of the characteristic vortices appear. 



applied in other phase field models such as the Cahn-Hilliard equation [TU] . 
Initial efforts to apply Newton-Krylov to the Ginzburg-Landau problem were 
taken in preconditioning is not discussed though. 

An important research topic in the context of the Ginzburg-Landau equations 
is the formation of vortex patterns in the solutions (see figures [T] [6cJ 6d ) . 
To understand the formation and dynamics of those patterns, the tools of 
nonlinear systems analysis can be employed. For example, numerical continuation 
techniques help computing a family of solutions as a function of a problem 
parameter, e.g., the strength of the externally applied magnetic field or the 
electric current at one of the boundaries. The main application of numerical 
parameter continuation is the construction of a bifurcation diagram that identifies 
the stability regions and the transition between stable and unstable patterns 
marked by bifurcation points |12j . A systematic bifurcation analysis of the 
patterns that appear in mesoscopic superconductors is carried out for square- 
shaped domains in [13] . The main computational load in numerical continuation 
are the linear solves with the Jacobian operator. By the sheer number of 
unknowns, this is particularly expensive for discretizations of three-dimensional 
domains. It is thus required to develop linear solvers for which the memory 
requirements and the computational cost grows slowly with the number of 
unknowns. To the knowledge of the authors, no linear scalable method for the 
Ginzburg-Landau problem has been developed. It is the goal of this paper 
to display that an AMG-preconditioned Newton-Krylov method is a viable 
approach for the extreme-type-II Ginzburg-Landau equations. 

The remainder of the paper is organized as follows. Section [2] reviews the 
Ginzburg-Landau equations for extreme-type-II superconductors; section |2.1 
is concerned with its linearization, the Jacobian, and discusses properties with 
respect to numerical algorithms. While section[3]introduces the applied discretiza- 
tion and shows that many important properties carry over from the continuous 
framework, section [4] is concerned with the solution of the Jacobian system and 
introduces a multigrid strategy. The convergence behavior is explored through 
numerical experiments on representative two- and three-dimensional domains. 
The document concludes with a discussion of the obtained results. 



2 



2. The Gibbs energy and the continuous Ginzburg Landau problem 



For an open, bounded domain Q C R 3 with a piecewise smooth boundary 
dfl, the Ginzburg-Landau problem is usually stated as a minimization problem 
of the Gibbs energy functional 



G(^,A)-G n = £^ f 
P Jn 



IVf + ^h/f + ll-iW-^f 



n 2 (V x A) 2 - 2k 2 (V x A) ■ H 



(1) 



over ip £ H^,(£l) and A <E ff^„(f2) [6 . The scalar- valued function ip is commonly 
referred to as order parameter, A is the magnetic vector potential corresponding 
to the total magnetic field. The physical observables associated with the state 
(ip, A) are the density pc — \ip\ 2 of the superconducting charge carriers (Cooper 
pairs) and the magnetic field B — V x A. The constant G n represents the 
energy associated with the entirely normal (non-superconducting) state. 

The energy ([I]) is presented in its dimensionless form, and it depends upon 
the impinging magnetic field Hq and the material parameters a, /3, A, £ G HL The 
ratio k := A/£ of the penetration depth A (the length scale at which the magnetic 
field penetrates the sample) and the coherence length £ (the characteristic spatial 
scale of ip) determines the type of the superconductor: It is said to be of type 
I if k < 1/^/2, and of type II otherwise. The two types behave fundamentally 
differently when exposed to a magnetic field: Type I superconductors exhibit 
alternating superconducting and nonsuperconducting regions, while type II 
superconductors show vortex patterns [14] (see figure [lj . 

Starting from the Gibbs energy and using standard calculus of variations, 
it is possible to derive the Ginzburg-Landau equations [TS] , a boundary- value 
problem in the unknowns ip and A. As anticipated in the introduction, we 
will simplify the problem and consider only the limit k — > oo (extreme type-II 
superconductors): this approximation gives satisfactory results for all high- 
temperature superconductors which have large values of k (typically 50 < k). 
In this case, the Ginzburg-Landau equations decouple for ip and A, such that 
the magnetic vector potential A is given up to gauging by the applied magnetic 
field H through 

VxA = H inR 3 , 

and V e X, X C Fg(fi), by 

f(-iV-A) 2 V-V>(l-|^| 2 ) infl, 
= GC(^) := { (2) 
[n-(-iV-A)V> on dn, 

with n being the outer normal on dVl and X :={?/> € Hc(Q) : GC(ip, A) bounded}, 
i.e., the natural energy space of Q. 

As the domain is given in units of £ = A/k, the large-K limit implies A ^> £ 
which means that Hq is not disturbed by the magnetic field induced by the 
electric charge density pc- 

Note that, for any given x€R, 

G£(exp(ix)ip) = exp(ix)GC(tp). 
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Consequently, any given solution ip of the Ginzburg-Landau problem is really 
just a representative of a whole set of solutions [ip] := {exp(ix)i/>: X € K}. 
This expresses the fact that for superconducting states, the actual value of the 
argument of ip is of no physical relevance: \ip\ 2 represents the observable. As 
the complex argument of any coefficient does not play any role in the scalar 
multiplication, c[ip] = \c\[tp], c G C, tp G X , it is natural to restrict the scalar 
field to M. The inner product in the vector space X over the field M is 



(0,^) R :=R<0,^)=fc / # . (3) 



n 



2.1. The Jacobian, the kinetic energy operator, and their properties 

Equation (J2| is a nonlinear equation in ip and hence classically suited for 
treatment with Newton's method. While there were efforts to solve ^ with a 
modified algorithm |16j , the generic approach of the full Newton system is applied 
here for its attractive second-order convergence. In this section, properties of 
the (continuous) Jacobian system 

J(iP)SiP = -&W>) (4) 

with 

J(il>)tp := ((-iV - A) 2 - 1 + 2|^| 2 ) ip + <p 2 7p. (5) 

will be discussed. Note that J(ip) is only linear if X is defined as vector space 
over the field R. 

The kinetic energy operator. Before analyzing the Jacobian operator J(ip) as a 
whole, we will take a close look at the part that is commonly referred to as the 
kinetic energy operator, 

Kip := (-iV - A)V (6) 

This operator is linear in X and self-adjoint with respect to the ordinary L 2 (f2)- 
inner product (see [T3]). Consequently, all eigenvalues of K are real- valued. Even 
more can be stated about its spectrum: From 

/ -0(-iV - A) 2 (p = [ (-iV - A)tp(-iV - A)<p - i / ^n-(-iV-A)^ 
Jn Jn Jon 

for all i/j,tp £ Lj.(f2), it follows that the kinetic energy operator is positive- 
semidefinite over the subspace X C X, 

I:={^el:n- (-iV - A)tp = a.c. on dfl}. 

This is because for all ip E X, 

(ij,K^) L2m = [ ^)Ki)<m= [ ||(-iv- a)^|| 2 dn>o. 

Jn Jn 
Moreover, the value of is attained if and only if 

(— iV — A)ip — a.e. on Cl, 
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from which in turn follows that 



= V x (-iVV>) - V x (Aip) = -i(V x W) - (V x A)ip - (Vip) x A 

= -Bi/j - i(Aip) x A = -Bip a.e. on ft. (7) 

Hence, only for vanishing magnetic fields B, the kinetic energy operator K is 
actually degenerate. 

An approximation for the smallest magnitude eigenvalue around the constant 
zero-field A = can be obtained by eigenvalue perturbation. Note that K(A ) 
is the Laplace operator with homogeneous Neumann boundary conditions, so the 
smallest magnitude eigenvalue of K(Aq) is Ao = 0, the corresponding constant 
eigenfunction v = 1. For the perturbed problem (K(A ) + SK)(v Q + 5v) — 
(A + SX)(vo + Sv), one gets 

6\{v , v ) = (v , (SK)vq) + (v , (5K - ISXjSv) , 

such that, in first-order approximation, 

sx ^ (vq, (SK)v ) = (v , (K(A) - K(A ))v ) 

Noting that K(A)v = (— iV - A) 2 v a = A 2 v 0l this yields 

A-Ao + ^^kW^. (8) 

This shows a lot more of the structure of the Jacobian operator J(i[>) al- 
ready: For any given ip £ X, J(ip) is the composition of a self-adjoint, positive- 
(semi) definite operator and some reaction terms. 

It is possible to infer certain properties of J starting from here. From a 
numerical point of view, insight into the adjointness and the spectrum of the 
operator will be highly desirable. The peculiar structure of J(ip), acting on ip 
and its pointwise complex conjugate, together with the inner product ^ in X, 
yield 

Lemma 1. For any given ip £ L 2 (f2), the Jacobian operator J(ip) is linear 
and self-adjoint with respect to the inner product ([3]). 

Proof. See Q3]. 

Now that the spectrum of J(ip) is known to be a subset of K as well, the 
natural question to ask is whether or not J(tp) is generally definite. Unfortunately, 
no such thing is true. Quite the contrary: Note that, for any solution ip s of j2j, 
we have 

J(ip s )(iip s ) = [(-iV - A) 2 - 1 + 2|V^ S | 2 ] (r0 s ) - 

= (1 - \A\ 2 ) (ty B ) ~ + 2i^ s 2 - i^Ws = 0, (9) 

and hence span{i?/> s } C ker J(ip s ). This is a direct consequence of the fact that 
QC(ip) (jij is invariant under the transformation ip = exp(ix)?/> for any x G M. 

Besides the fact that there is always a degenerate eigenvalue and that all 
eigenvalues are real, not much more can be said about the spectrum; in general, 
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J{tl>) is indefinite. The definiteness depends entirely on the state tp; if if) is a 
solution to ([2]), it is said to be physically stable or unstable depending on whether 
or not t7(V0 has positive eigenvalues. Typically, solutions with relatively low 
energy tend to be stable whereas solutions with relatively high energy tend to 
be unstable. 



3. Discretization in finite volumes and link variables 

In recent years, the research in applications for superconductors has taken 
strong interest in studying the effect of the sample geometry on superconductivity 
phenomena, for example, of dents or holes in a domain. Such geometries cannot 
be captured well by classical Cartesian staggered grids [13], so finite element 
and finite volume approaches have been developed that incorporate properties 
of the continuous Ginzburg-Landau equations such as the gauge invariance. In 
|17j . the method has been described for two-dimensional domains and shall be 
described here in general terms. 

Let Xj £ R d , d £ {2, 3}, j £ {1, . . . , n} be a given set of discretization points 

at which states tp £ X will be approximated by ipi h ' ~ 4>{ x j)i G C™. Each 
discretization point Xj be equipped with its corresponding Voronoi region, 

Vj := {x £ M. d : \\x - Xj\\ < \\x - x k \\ Vk ^ j}. 

The set {V}}™ =1 is referred to as Voronoi tessellation corresponding to the 
generator set {xj}" =1 . The dual to a Voronoi tessellation consists of simplices 
and is referred to as Delaunay triangulation {T i \™ =1 (see, e.g., figure[3]). 

For the domain Q,^> := (Jj=i = U£li the significant part of the Gibbs 
energy ([!]) can be written as 



F(^,A) = ^/ ||-iV^-A^|| 2 + ^ / -|^| 2 + 2 

i=l jTi 7=1 Vj ^ 



2 H 4 



The second term, F%, is readily discretized by mass lumping, 

*i h \* w ) ■■= E m (M h) \ 2 + livfi 4 ) ■ cm) 
3=1 ^ * 

For the discretization of F\ , we will first refer to a technique for triangular 
meshes in |17J . extended to arbitrary dimension here. 

Lemma 2. Let e^. i £ {1, . . . ,n}, with n := d{d + l)/2 be the edges of a 
nondegenerate d- dimensional simplex. Then the symmetric rank-1 matrices 
{eiej}\ l =1 form a basis of the vector space of symmetric d x d-matrices. 

Proof. The number n of edges in a d-dimensional simplex coincides with the 
dimensionality of the vector space of symmetric d x d-matrices. Hence, only 
linear independence has to be shown. Assume then that 

n 

o d , d = J2^( e i e I) (u) 
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with some j3 € R™. Since the simplex is not degenerate, there is a regular matrix 
U that maps the edges {ei}™ =1 onto the edges {ei}™ =1 of the unit simplex. With 
this, (11) is equivalent to 

(n \ n n 

Y^^eJ) U T = Y,fc(Ue l (Ue l ) T ) =^00^). (12) 
i=l / i=l i=l 

For the edges parallel to one of the axes, we have e-iej = ete^ with being the 
unit vector in /c-direction. For the edges between the two axes k%, ki, we have 
e-iej = (efej — e-k 2 )i e k 1 — ek 2 ) T ■ As the matrix that belongs to edge between fci, 
ki is the only matrix with a nonzero entry at (fci, k-i) (namely —1), its coefficient 
in ( 12 ) must be 0. Similarly, the same holds for all other coefficients, such that 
( [II I can only be fulfilled of /3j = for all i £ {1, . . . ,n}. Hence, the matrices 
{eiej}2 =1 are linearly independent. 

Since the {e i e^ r }™ =1 form a basis, there exists in particular a unique set of 
coefficients a 6 K n , such that 

Id,d = E a i( e i e I)- 

edges Gi 

From this, we immediately conclude 

Corollary 3. Given a nondegenerate simplex S € K d with edges e^j :— Xi — Xj, 
i,j G {1, . . . , d + 1}, i 7^ j, there are coefficients ai_j such that 



j \\u\\l = \S\ (U,U) 2 = ^ a hi < U ' e M>2 ( e i,j> U )2 = E a ij\( e i,3i U ) 



edges eij edges 

for any u E C d . 

Given a simplex, one way of determining the edge coefficients a t j is to solve the 
symmetric and positive-definite linear equation system Ma = b with 

Mij := {ei,ej)l , h := |5| He,]^ , 

where the edges are indexed subsequently. 

Remark 1. For triangles, the edge coefficients ctij are explicitly given by 

1 ,a 1 **J 

CXi j = - COt (fi j = , 

J 2 J 2 

1 



ej.fc e,- k 



where 9i ,• is the angle opposing the edge e< ,• [T7], and tj , ■ := , 
with fc ^ {i, j}. 

With corollary [3] (and the coefficients a, j from there), i*\ can be approximated 

by 

m 

F 1 (i;,A)^F 1 (if,A):=J2 E «S|e^-HV-A)^(% fc )| 2 

2—1 cdgCS Bj fc of Tj 
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with Xj^ := \{xj + ccfc), or, more compactly, 



F 1 (tp,A)= ^2 a iA e 3,k ■ (~iV - A)ip(x J:k )\ 2 



(13) 



edges e 7 - 



with the edge coefficients 



a j,k 



E 



simpliccs Ti 
adjacent to edge ej k 



One could now do a finite difference approximation in the differential terms 
of ( 13 ) to receive a Gibbs energy defined over the discretized function space 
C". Note, however, that such naive discretization schemes of the momentum 
operator — iV — A lead to systems that preserve gauge invariance - inherent 
to the Ginzburg-Landau equations - only up to a certain order in the spatial 
discretization. It is hence customary to rewrite the momentum operator in 
terms of variables that ensure preservation of gauge invariance for any pointwise 
discretization. Following |16) . for any given normalized spatial direction v, let 



U v (x) := exp ^— i J v ■ A(w) dw 



(14) 



with arbitrary, fixed x E x + sp&n{v} (e.g., x = x — (x ■ v)v). Since U v (x) sits 
on the unit circle, one has U v (x) U v (x) = 1, and with this 



U v v ■ V(U v i>) ee U v (-iv ■ AU v ip + U v v • Vip) = iv • (-iV - A)ip. 
Thus, F\ can be written as 



(15) 



edges e Jifc 



U ejik (x j>k )e j>k ■ V(U ej k ip)(x^ k ) 



Finite difference approximation finally yields the discretization 



with 



edges e jt k 
edges e 3 - jj 



Uj t k '■= exp — i / ej k ■ A(w) dw 



often called link-variable |16j . If A is known only at certain points of along the 
edges, Uj t k could again be approximated by a quadrature formula. 

Finally, together with (10), the discrete Ginzburg-Landau energy functional 
is defined as 

FW^^U) := F^ h \^ h \A) + F^(^). (16) 
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The standard Euler-Lagrange formalism now yields a necessary condition for 
extremal points of the energy functional, 



= 23?| 



E a ^ k 

, edges e j:k 



This is equivalent to the discrete Ginzburg-Landau equations, 

V* G {l,...,n} : = (i^W) . - (l - |^ (/l) | 2 ) , 

where the discrete kinetic energy operator is defined by 

[4 h) - u ^i h) ) 4 h) + - 

with the discrete inner product 

n 



E a ^ 

edges e Jifc 



,.('0 



(17) 



(18) 



(19) 



Remark 2. In matrix form, the operator K^ 1 ' is represented as a product 
KW = D~ X K of the diagonal matrix D^ 1 , Z?^ = \Vi\, and a Hermitian matrix 
K. 



The discretization (17) has several advantages, starting with the fact that 
the boundary conditions of the Ginzburg-Landau equations ^ are naturally 
contained. Also note that the discrete kinetic energy operator ( |18[ ) coincides, 
up to the terms Uj t k, with the discretization of the Laplace operator with 
homogeneous Neumann boundary conditions. Similarly, it has a number of 
desirable properties that will make the iterative solution of the Jacobian system 
easier. 



Lemma 4. The discrete kinetic energy operator (18) is self-adjoint with 
respect to the discrete inner product (19). 

Proof. Let <f>W,^ h ) e C™. Then 



jfCtyCO 



E a ^ k 

edges e jtk 



(+5 



E 



atj,k 



{4 



Lemma 5. The discrete kinetic energy operator (18) is positive-semide finite. 
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Proof. Let € C". Then 

edges e i k 

= E 

edges e 3ifc 



Noting that Uj^Uj^ = 1, this yields 

>,W ir(h),i,(h)\ 



ay,* 



> o. 



(20) 



For Aq = 0, the state -tp^ = 1 is obviously an eigenvector of IfW with 
the eigenvalue 0. Equation (20) also delivers an approximation Ao for the 
smallest-magnitude eigenvalue for perturbations of Aq, namely 



A = 



edges e jik 



in^r 1 E «a* 

edges e 3jfc 



2 sin 



or, in first approximation, 



edges e 3ifc 



H^r 1 E 

edges e 3ifc 



e jt k ■ A(w) dw 



(21) 



Compare this with the corresponding continuous expression ([8]). 

Completely analogous to the results for the continuous Jacobian operator 
J(ip), the discrete Jacobian operator 

(jWtyW)4>W) := (#(*V*>) + (-1 + 2|^f>| 2 ) + (vfV^ 



of (17) is self-adjoint with respect to the inner product 



(22) 



Like the continuous Jacobian operator J(V0, J^ h '(if)^) also has a nontrivial 
kernel if V'*'' 1 ' 1 is a solution to the problem. While in the Newton process, the 
Jacobian system will never need to be solved in exactly a solution, states close 
to a solution might introduce numerical difficulties when nearly-singular systems 
need to be solved. Techniques for this situation include adding phase conditions 
[T3l or deflation methods. 
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Remark 3. Note that there is a vector space isometry of C n as vector space 
over the field K and R 2 ™ with its natural inner product: For all <f>W , ip^ £ C", 



This is relevant in practice if the complex- valued original problem (17) in C™ is 
implemented in terms of M. 2n . Using the natural inner product in this space will 
yield the expected results without having to take particular care of (v) R . 

4. Algorithm and numerical results 

For real-world three-dimensional domains, the solution of the discrete equiv- 
alent of the Jacobian system Q will have too many unknowns for black-box 
strategies such as LU -decompositions to fit into memory. Exploiting the sparsity 
structure of the operator is crucial, and hence Krylov subspace methods are an 
attractive alternative. The choice of the appropriate Krylov subspace method 
depends on the spectral properties of the operator. Its self-adjointness (see 
lemma [l] and its discrete equivalent) make it possible for symmetric Krylov 
subspace methods to efficiently solve the linear system if adapted for the inner 
product ( |22[ ) (see also remark [3]). This avoids the larger memory- requirements 
that come with asymmetric solvers such as GMRES. Furthermore, as J^ h \^ h ^) 
is generally indefinite (depending on and the number of negative eigenvalues 
can be large, CG may be unsuitable as a solver. While it is known to perform 
well for indefinite problems when the number of negative eigenvalues is not too 
large [18] . convergence can be irregular. In contrast, MINRES is designed to 
deal with indefinite systems and is hence a more suitable choice. 

4-1. Preconditioning 

As the main computational effort of the nonlinear solver flows into the linear 
solves of the Jacobian system, and the complexity of the linear solve usually 
grows faster than linearly with the number of unknowns in the system, it is 
crucial to explore the possibilities of accelerating the Krylov solver using an 



appropriate preconditioner. Given the results of section 2.1 we will evaluate the 
use of approximate inverses of the operator 

as a preconditioner for JW^W). The operator PW^C 1 ') is obviously self- 



adjoint with respect to the standard discrete inner product (19) and positive- 
semidefinite. From lemma [5j we can conclude that it is even strictly positive- 
definite except for the uninteresting case ip^ = 0, A = 0. Moreover, 
is derived from a geometric discretization, and its sparsity structure coincides 
with that of the Laplacian with homogeneous Neumann boundary conditions. 
This makes the inversion of P^ h \tp^) a suitable target for algebraic- multigrid 
(AMG) strategies which are known to yield optimal convergence behavior in the 
sense that the number of iterations required to reach a certain stopping criterion 
is independent of the number of unknowns in the system. Furthermore, AMG 
methods are memory-efficient and scale well in parallel computing environments 
[Hfl |2T)1 |2"T] . The only nonstandard circumstance here is the fact that the matrix 
entries are complex-valued. Difficulties in this area, however, were discussed and 
treated in |2"2"1. 
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Remark 4. The operator Q( h \ip( h '>) defined by 

= (kW + 2|?A ( ' l) | 2 ) <p (h) + {4> {h) ) 2 W ] V<t> {h) G C" 
is obviously self-adjoint and also positive-definite since 

[^ h \2\^\ 2 < t>^ + {^) 2 4> {h) ) 

n 



i=l 



En 



i=l 
> 0. 



/ CO j^) 



It would hence also be a candidate for a good preconditioner. However, unlike 
P^ h \ it cannot be represented as a matrix and is thus not suitable for solution 
with AMG. 

Note that the operator AMGfc (A, b) , defined by k AMG cycles applied to a 
Hermitian problem Ax = b, is again Hermitian. With remark [2] (and D from 
there), we have that 



(hy 



= (k + 2D\^\ 2 ) D, 



so the approximate inverse of pW(ipW), 

R [ k\^ {h) ) = AMG fe (k + 2D|V> 



(ft) 1 2 



D- 



(23) 



is self-adjoint with respect to the standard discrete inner product (19). 

We will now explore this idea through numerical experiments, for the precon- 
ditioners R± {ip^ h%> ) and 



R&H^) '■= (P^^)) -1 = (k + 2D\^\ 



D, 



(24) 



where K + 2D\\p^ h '\ 2 is inverted numerically with high accuracy. 

It is notoriously difficult to rigorously characterize the spectrum of the Ja- 
cobian operator of the Ginzburg-Landau problem, and the situation is sim- 
ilar for the preconditioned operator. Nevertheless, if (A, <^>C0) is an eigen- 
value/eigenvector pair of the preconditioned operator Roo(ip^) jCO^CO), i.e., 

one gets 



A = 



(0W,PW(^W)^{*)) 



= 1 



3?((^) 2 ,(^W) 2 



(25) 
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10 -14 

10 20 30 200 400 

number of MINRES iterations relative complexity 

Figure 2: Typical residual behavior for MINRES, applied to the problem 
jWtyW^h) = b (h) f preconditioned with (tp ih) ) for different k, with each 
ip( h ) = 5C 1 ) = i ) initial guess (j)^ = 0, for fli'q in A z . The number of unknowns 
is 1000 2 in all cases. 10 AMG steps solve the preconditioning problem in each 
MINRES step up to at least ||r p || < 10~ 12 . In the figure on the right, the 
computational complexity is measured in terms of the cost of one matrix-vector 
multiplication. For this setting, the application of one V-cycle costs as much as 
about 3.31 matrix-vector multiplications. 



In case > 1 (which can happen during the Newton iteration), the eigen- 

values cluster around 1 ± \ (depending on the sign of ?ft((4>^) 2 , (tp^) 2 )), so 
the preconditioned problem can be expected to be solved in a small number of 
Krylov iterations. Noting that solutions tp of the Ginzburg-Landau equations 
([2]) fulfill \tp\ < 1 pointwise, though, ( |25| ) unfortunately gives little insight in 
the behavior close to a solution. The same is true for the bounds gained from 
estimating the denominator term {(jr- h \ K^ h ' <p^ h n with the help of the smallest 
eigenvalue approximation for weak fields ( |21| . 

While R&\ipW) is obviously more expensive to apply, it is expected that it 
will yield a smaller number of Krylov iterations as compared to preconditioning 
with R[ h) (il)(V). Figure |] illustrates this: For a fixed setup, the preconditioners 
R^\ip^) with k G {1,2,5,10} are compared, where for this particular case 

^Hoi^^) ~ Rooti^) t° machine-precision. Preconditioning with R^q (ipW) 
indeed results in the smallest number of required MINRES iterations; if fewer 
V-cycles are applied per iteration, the number of iterations increases. A better 
measure for the overall computational cost than the sheer number of Krylov 
iterations, however, is the number of performed V-cycles together with the 
matrix-vector products. While the latter mainly depends on the number of 
nonzeros in the kinetic energy operator K^ 1 ' , the cost of the former also depends 
the many parameters of AMG. In all of the experiments performed in this paper, 
the cost of the application of one V-cycle is between three and four times the cost 
of a matrix-vector product of the corresponding matrix. As can be seen in the 
right panel of figure [2j no more than the equivalent of about 140 matrix- vector 
products is are required in total to converge the MINRES process in combination 
a single V-cycle preconditioning. At the same time, 10 cycles per step require 
the equivalent of about 480 matrix-vector multiplications. This points to the 
fact that the approximate inversion with a single V-cycle will lead to the fastest 
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solver. 

For the experiments in figure [2] and all experiments in the remainder of 
this paper, smoothed-aggregation AMG with one pre- and one post-smoothing 
step of symmetric Gauss-Seidel was used. The method is implemented using 
PyAMG [23]. 

We now look at the application to the two-dimensional regular polygons in 
the 2-y-plane 



-5^3/2 
-5/2 



"tri • -^convc 

il sq := {x : \\x\\ x < 5/V2} (figure 
^circ := {x : \\x\\ 2 < 5} (figureptel 



5^/2 
-5/2 



(figure 3a) 



3b I 



and the three-dimensional regular polyhedra 



9. 





(figure 3d), 



2^/3 


-1/3 

^cubc := {x : WxW^ < 5/V3} (figure[3| 
HbaU : = {x ■ \\x\\ 2 < 5} (figure |3§, 

all centered at the origin with circumradius 5. For each domain, both the 
potentials 

A 9 (x):=±(-y,x,0) T , (26) 
representing the homogeneous field B = (0, 0, 1) T , and 

A d (x) := pz^ipr On x (x - x )), 

representing the inhomogeneous field generated by a magnetic dipole at the 
location 

' (0, 0, 1) T for the 2D domains, 
(0,0, 6) T for the 3D domains, 



x Q 



and with the dipole moment m — (0, 0, 1) T , are considered. For all experiments, 
we considered JW(^W)) with = 1. For other choices of see the 
paragraph on numerical continuation below. 

Figures [2] and [5] show the number of MINRES iterations as a function of the 
dimension of the solution space. For the unpreconditioned system, the number of 
iteration increases as expected since the finer discretization makes the condition 
number of and hence larger. In contrast to this, when and 
are applied as preconditioners, the number of iterations remains bounded in all 
problem settings as the discretization refines. Although the number of iterations, 
when preconditioned with , is slightly larger compared to preconditioning 

with R^ , the former is actually computationally cheaper as discussed above (see 
figure [2]) . These numerical experiments suggest that for various fixed domains 
and magnetic vector potentials, the number of iterations of the Krylov solver is 
independent of the number of unknowns. 
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Numerical parameter continuation. A common application of the Newton-Krylov 
solver is in a numerical continuation context where a family (or curve) of solution 
states (from a given state space X) is constructed as a function of a parameter 
(from a given parameter space A) in the system. This is a popular way of 
exploring the solution landscape of nonlinear equations, and amongst the most 
widely used algorithms for this purpose is pseudo-arclength continuation j!2j . a 
predictor-corrector method. Here, in each continuation step an initial guess is 
constructed as an extrapolation to the solution curve in X x A which is then 
corrected perpendicularly to the extrapolation, typically involving a Newton- 
Krylov process. In applications, many curves, each with thousands of solutions, 
are computed. Each continuation step requires the solution of a nonlinear system, 
each of which requires the solution of a Jacobian system. 

As this setting presents a typical use case for the preconditioner introduced 
above, a representative problem is discussed in this section. As opposed to all 
previous numerical experiments, the state ip^ deviates significantly from the 
initial state tpQ = 1 in the corresponding numerical experiment (see figures 



6c 



6dj). 

Figure]^ illustrates the performance of the ^^■'-preconditioned Krylov-solver 
for f2sq' with edge length 10 in fiA z (26), [i G R. The strength [i of the magnetic 
field is used as continuation parameter, and the continuation is started with the 
trivial solution ip^ = 1 at /x = 0. As /i increases, the solution starts to deviate 
from the homogeneous state. Throughout the parameter continuation, vortices 
appear in the domain and the state loses its stability [13] . a process marked 
by eigenvalues of the Jacobian crossing the origin, i.e., a change of definiteness 
of the Jacobian operator. The right panel of figure [6] shows, for each point on 
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Figure 4: MINRES performance for different two-dimensional domains and 
magnetic vector potentials at the state i/jq = 1. The plots show the number of 
iterations necessary to reach the relative residual of 10~ n in the norm given 



by ( 22 1 as a function of the dimension of the problem. Starting guess for the linear 



iterations is_^p = throughout. Plotted are re sults for the unpreconditioncd 



problem 
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[ :. and the preconditioner 
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Figure 5: MINRES performance for different three-dimensional domains and 
magnetic vector potentials at the state i/jq = 1. The plots show the number of 
iterations necessary to reach the relative residual of 10~ n in the norm given 



by ( 22 1 as a function of the dimension of the problem. Starting guess for the linear 



iterations is 4>q = throughout. Plotted are re sults for the unpreconditioncd 



problem 
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the continuation curve, the number of iterations in the preconditioned MINRES 
solver that was required to solve the Jacobian system up to ||r|| < 10~ 8 in the last 
Newton step. While initially around 50 iterations are required, the introduction 
of an unstable eigenvalue at the swallow tail (fi s» 0.30) slows down MINRES 
convergence. This is due to the fact that the positive-definite preconditioner 
does not capture negative eigenvalues. While the convergence is not slowed 
down by an order of magnitude, the highly unstable high-energy states around 
fi = 0.8 require up to 150 MINRES iterations to converge. The local peaks in 
the iteration requirements are due to the inherent loss of orthogonality of Krylov 
basis vectors in MINRES [21] ■ This effect can be alleviated by storing and fully 
reorthogonalizing the Krylov basis in each MINRES step such as implemented 
in GMRES. 



5. Conclusions 

The authors propose a preconditioned Newton-Krylov algorithm that solves 
the extreme type-II Ginzburg-Landau equation. The solution method uses an 
AMG preconditioning strategy that yields optimal convergence and scalability 
for mesoscopic domains. 

The Ginzburg-Landau operator consists of a kinetic energy operator that 
depends on the given magnetic field and a nonlinear reaction term. The lineariza- 
tion of the operator is analyzed and it is found that the Jacobian of the system 
is self-adjoint with respect to the nonstandard inner product (J3j) . Its spectrum 
is indefinite if ip describes a physically unstable solution of the equation. The 
properties of the kinetic energy operator K are also discussed and it is found 
to be self-adjoint and positive (semi)-definite. These properties are maintained 
after discretization with finite volumes and link variables. The proposed precon- 
ditioner takes advantage of this by applying an algebraic multigrid scheme to 
the operator pW^W) = + 2\tp^\ 2 . Numerical results for representative 
domains point towards the optimality of the algorithm in the sense of indepen- 
dence of the number of linear solver iterations from the discretization resolution. 
This suggests that, qualitatively, no further improvement can be reached. 

Moreover, the performance of the preconditioner is assessed in a numerical 
parameter continuation context where a family of solutions is generated for 
changing strength of the applied magnetic field. The good convergence results 
from the test domains are confirmed here. The presence of negative eigenvalues, 



however, slows down the Krylov convergence if used with a CG solver (figure 6b). 
Moreover, other factors, such as a large domain size, have shown to hamper 
the convergence. To gain deeper insight into the convergence behaviors, clearer 



results than (251 on the spectrum of (pW) _1 jW are needed. 

Nevertheless, this research opens up new possibilities for the exploration 
of the energy landscape of type-II superconductors. Computation of three- 
dimensional problems are now accessible with grid resolutions on par with 
current two-dimensional calculations. 

A natural extension of the presented work is to approach the solution of the 
full Ginzburg-Landau problem in which the magnetic vector potential cannot 
be treated as given [15j. Numerous numerical and computational challenges are 
posed there, e.g., how to efficiently solve the Jacobian system. The presented 
preconditioner could be used to construct a block-preconditioning strategy for 
the general (nonextreme-type-II) Ginzburg-Landau equations. 
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(a) Gibbs energy as a function of /i. Note 
the so-called swallow tails which are typ- 
ical for numerical parameter continuation. 
Structural transitions in the solutions usu- 
ally happen around here. 




(b) Number of MINRES iterations till ||r|| < 
10 — 8 in the last Newton step as a function 
of fi. With increasing energy, stability is 
lost and the presence of unstable eigenmodes 
hampers the performance of the precondi- 
tioner. 





(c) Solution \ip\ 2 , argV at continuation step 22 with fi a 0.47, w -0.47. The 

system has undergone a first transition and four vortices have moved in. 




1 




o 

(d) Solution |»/>| 2 , arg^ at continuation step 49 with fj. sa 0.93, F^ « -0.12. The 
system has undergone a second transition and now contains eight vortices. 

Figure 6: The performance of the preconditioned Krylov solver in the context of 
numerical continuation in [iA z for f2sq with edge length 10, 1000 2 unknowns. 
The solution is continued in the parameter [i with the help of pseudo arc-length 
continuation where each continuation step requires the solution of a nonlinear 
system. 
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